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Current status of Dynamical Overlap project 

N. Cundy a 

a Theoretische Physik, Universitat Wuppertal, Gaussstrasse 19, D42109 Wuppertal 

We discuss the adaptation of the Hybrid Monte Carlo algorithm to overlap fermions. We derive a method 
which can be used to account for the delta function in the fermionic force caused by the differential of the sign 
function. We discuss the algoritmic difficulties that have been overcome, and mention those that still need to be 
solved. 
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1. INTRODUCTION 

The overlap operator pQ is the closest known 
lattice Dirac operator to the continuum operator. 
With the lattice community currently moving be- 
yond the quenched approximation, one obvious 
possibility is to use dynamical overlap fermions. 
The advantages of the overlap operator are well 
known: it satisfies the Ginspag- Wilson lattice 
chiral symmetry exactly; there is an easy non- 
perturbative renormalisation; there is no opera- 
tor mixing involving different chiral sectors; there 
is a well defined index (Qf — ^Tre(Q)), equal 
to the topological charge in the continuum limit; 
the anomaly is correctly accounted for; and it is 
essential for stuties of topics such as topology, 
spontaneous chiral symmetry breaking, and the 
eigenvalue spectrum of the Dirac operator. Since 
chiral symmetry is so important to low energy 
QCD, it seems a waste not to use the only Dirac 
operator which we know of which fully respects 
this symmetry 

Of course, there are reasons not to use the over- 
lap operator. Firstly, it is considerably slower 
than (for example) staggered or clover fermions. 
Secondly, the discontinuity in the overlap opera- 
tor creates a number of unique problems when 
trying to impliment a Hybrid Monte Carlo al- 
gorithm. The first of these issues will not be a 
problem once we have sufficiently fast computers, 
which will be in the very near future (we have 
already run some trajectories on a 16 3 32 lattice). 
Now is the perfect time to tackle the second prob- 
lem, and to create a Hybrid Monte Carlo algo- 



rithm for overlap fermions. 

In this talk, I shall summerize the work done 
so far on this issue. Work has been published 
in this area by Z. Fodor et al |2I3I4| . by myself 
in collaboration with Thomas Lippert and Stefan 
Krieg [EEEI, and by T. DeGrand and S. Schae- 

fer iHEma. 

Section |21 provides a brief introduction to Hy- 
brid Monte Carlo and the overlap operator. Sec- 
tion[3]discusses the problem of topological charge 
changes. Section 0] mentions some additional 
problems and advantages concerned with dynam- 
ical overlap fermions. Section gives a few nu- 
merical results, and our conclusions are presented 
in section 

2. HYBRID MONTE CARLO WITH 
THE OVERLAP OPERATOR 

The overlap Dirac operator is 

£> = (l + /i)+76(l-/i)e(Q), (1) 

where /i is a mass parameter, Q is the hermitian 
Wilson Dirac operator, and e is the matrix sign 
function. In our numerical simulations [H] , we use 
a Zolotarev rational approximation to the sign 
function with the small eigenvalues treated 
exactly using eigenvalue projection, but for the 
purposes of this talk, I shall assume that we can 
calculate all the eigenvalues and eigenvectors 
{ipi | of Q and thus treat the matrix sign function 
exactly using e(Q) — sign(Ai) This will 

simplify the algebra while retaining the important 
features of the algorithm. Wc will define H as the 
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Hermitian overlap operator 75-D. 

The Hybrid Monte-Carlo (HMC) algo- 
rithm updates the gauge field in two steps: 
(1) a molecular dynamics (MD) evolution of the 
gauge field; (2) a Metropolis step which renders 
the algorithm exact. In the MD step, we in- 
troduce a momentum, II, which is conjugate to 
the gauge fields U, and a spinor field which 
is used to estimate the fermion determinant via 
a heat-bath. We define an Hybrid Monte-Carlo 
energy 

E= 1 -Il 2 + S g [U]+<p\H)- 2 <p. (2) 

S g is the gauge action, and we will use either the 
Luscher-Weisz or Wilson plaquette action. We 
introduce a computer time r and integrate over 
the classical equations of motion to generate the 
correct ensemble. We cannot perform an exact in- 
tegration, so we need to use a numerical method 
such as the Omelyan integration step This 
will create a small error in the energy conserva- 
tion, which we can correct for by including an 
additional metropolis step, accepting or rejecting 
the new configuration according to a probability 
P acc = min(l,exp(A), where A = Ei - E f , E, 
is the initial energy and Ef the energy at the 
end of the MD. It is therefore important that 
we conserve energy as well as possible during the 
MD to ensure a high acceptance rate. Note that 
we do not have to use the classical trajectory: 
any reversible update which leads to a small A 
will suffice. Our MD procedure does not even 
have to conserve area, as long as we can eas- 
ily calculate the Jacobian J. If we have a non- 
area conserving update, we just have to include 
the Jacobian in the Metropolis step, using a new 
A = Ei — Ef + log J [T3|. To have a high accep- 
tance rate, we need e A ~ 1. 

The crucial part of the MD procedure is the 
force used to update the momentum, defined as 

_ _ rr 0(S g [E/] + <ftt(£)tj))-i0) 



The fermionic part of this force for the overlap 



operator is 

F F IL + ILFl = -(1 -11 2 ) 

^(*i<Q + w+ (4) 

We can differentiate the eigenvectors and eigen- 
values of Q using a procedure analogous to first 
order time independent perturbation theory in 
quantum mechanics jSJ. This gives 

- 

£ 1*) <*i § m in 'J^Mjz^M 

+ £> < ><^| JW(A0. (5) 

i 

Note that only mixings between eigenvalues of 
different signs contribute to the fermionic force, 
and only mixings between the small eigenvalues 
are important. The main feature of the fermionic 
force is the Dirac ^-function coming from the dif- 
ferential of the sign function. We shall discuss 
how to deal with this in the next section. 

3. EIGENVALUE CROSSINGS 

The (5-function in the fermionic force, which 
occurs whenever one of the Wilson eigenvalues 
crosses zero (i.e. whenever there is a change in 
the index of the overlap operator), should, in an 
exact integration, introduce a discontinuity in the 
momentum, which will exactly cancel the discon- 
tinuity in the pseudo-fermion energy caused by 
the abrupt change in the matrix sign function. 
We can visualise this in a classical mechanics pic- 
ture by picturing a potential wall of height —2d 
surrounding each topological sector. We can eas- 
ily calculate the height of the wall (which can be 
either positive or negative), either by integrating 
the fermionic force across the (5-function, or by 
calculating the difference in the pseudo-fermion 
energy. Both these procedures give 

d=-(l- f i 2 ) 

M(]^{75,e(A-)MM} (£^2 l*>- 

(6) 
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H + is the Hermitian overlap operator just af- 
ter the crossing, H~ the operator just before the 
crossing. In a classical mechanics picture, the mo- 
mentum would be updated in a direction parallel 
to 77, the normal to the topological sector wall, 
thus: (n+,?7) 2 = (n~,r7) 2 + id. If the momen- 
tum is too small (i.e. this procedure can lead to 
an imaginary momentum) we reflect of the po- 
tential wall. However, as remarked earlier, there 
is no reason why we have to stick to the classical 
mechanics picture. Below, we shall describe how 
a general updating procedure can be derived to 
account for the potential wall. 

To calculate the Jacobian, we need to work 
with the coordinate and momentum vectors u and 
7r. 7r can be calculated easily from the momen- 
tum field n M (a;) = T^ 1 (x), where u and x refer 
to the direction and lattice site, and Tj are the 
generators of the gauge group. The gauge co- 
ordinate u is defined so that an update of the 
gauge field U — > e m U corresponds to u — > u + n. 
We correct for the (5-function in three steps. (1) 
We update the gauge field to the potential wall 
u c = u~ + t c 7t~ ; (2) We update the momen- 
tum, using (7T+) 2 = (7rr) 2 + Gi(n~,u c ), where 
we shall determine the functional form of G later; 
(3) We return the gauge field to the original point 
u c = u + — t c 7t + . Here r c is the computer time at 
which the eigenvalue is zero, i.e. 



(u c - U,7]) 



(7) 



(tt~ ,rj) refers to the scalar product of the two 
vectors. Differentiating t c with respect to u and 
7r gives 



bian. Using the gauge update above gives 



J 



dut 



9tt+ 
duf 



dw+ 

a S 

duf 
du k 

dr c 
dir7 



'k 

duf dr c 



dul 



dul 



(10) 



We can calculate J using two quick determinant 
manipulations. We subtract r c times the top row 
from the bottom row. We then subtract t c times 
the right column from the left column. These 
manipulations kill the bottom left hand clement. 
Two lines of algebra later, and we obtain 



J 



(ii) 



We consider the components of the momentum 
normal to r\ and perpendicular to rj separately. 
For example, if we update the momentum normal 
to rj, then we use the momentum update above 
and insert this Jacobian into the condition e A = 1 
needed for a high acceptance rate. This immedi- 
ately gives us a differential equation for G v : 



-G v /2+2d 



1 



;t _ | 1 0G V \ 7T+ 



1. 



(12) 



It is trivial to solve this equation to obtain the 
momentum update 



-«) 2 /2 



-(7r-) 2 /2-2d 

A(|d|)(l-e 



-2d\ 



0. 



(13) 



8tc_ _ chc_ 

T fe c du k 



Vk 



dir k C du k c {it, r/) 



(8) 



Any function, g, of u c (such as d or 77) will obey 
the relation 



dg = dg_ 

duk C du k 



(9) 



We are now in a position to write down the Jaco- 



We have written the constant of integration as 
A(l — e~ 2d ) to ensure reversibility. A should lie 
in the range < A < 1. A = gives us the 
classical mechanics solution. 

We cannot allow the final momentum to be 
complex. Therefore, if the initial momentum is 
in the range a < exp(— (n~) 2 /2) < b, where 



„2d\ 



a=A(l- 



b=e zd {l-A)+A 1 



(14) 
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we can use equation (|13fl to update the momen- 
tum (we call this case transmission). Following 
we reflect the momentum of the potential wall 
if the momentum lies outside this range (i.e. we 
use ttJ = — 7r~, with some additional terms - 
outlined in - to ensure 0(r 2 ) energy conser- 
vation) . We need to keep the transmission rate as 
high as possible in order to reduce the topologi- 
cal autocorrelation. We can calculate the trans- 
mission rate by assuming that the momentum is 
initially distributed according to exp(— (n~ ) 2 /2). 
The probability of transmission is 



V"21og(b) 
\/-21og(a) 



(15) 



and we can show that dPt/dA > for > A > 1. 
Therefore A = 1 maximises P t . 

We can also update the momentum in direc- 
tions perpendicular to r/. The procedure is ex- 
actly the same as above, and there are a large 
number of possible solutions, depending on the 
number of dimensions and how we combine the 
gauge fields. The simplest is, perhaps, the two 
dimensional case. We write m = rcos9 and 
7T2 = rsin#. We can now proceed to change r 1 , 
using 

e -(r+) 2 /2 _ e -(r-) 2 /2-2d 

~A{\d\)(l-e^ d ) = 0. (16) 

We can also update as many of these pairs of mo- 
mentum fields as we like, or work in more than 
two dimensions. However, no matter what up- 
date we use perpendicular to the probability of 
transmission remains (min(l,e 2d ). However, we 
can use the updates perpendicular to r\ to remove 
one rather large annoyance: an 0(r c ) energy vi- 
olating term. Without some means of removing 
this, we would have to reduce the time step to 
unfeasibly small values: it would cripple the al- 
gorithm. 

It is easy to show that the procedure outlined 
above only conserves energy up to order r c , be- 
cause we update the gauge field, but not the mo- 
mentum field, to the topological sector wall. We 



can update the momentum field using the force 
perpendicular to n without affecting the Jacobian 
(as long as we only change components of the mo- 
mentum field orthogonal to this force), but we 
cannot so easily use the component of the force 
normal to n. This leaves us with an energy viola- 
tion AE = T c (F+,r))(r),Il+) - T c (F-,r))(r),Il-). 
We can, however, add this term to the perpendic- 
ular update, for example by changing d — > d—AE 
in equation Ijltil) . This allows us to remove the 
O(t) and many of the 0(r 2 ) energy violating 
terms, leaving us with a correction step that is 
almost 0(t 3 ). 

4. OTHER ISSUES 

We had 

rf = -(l-M 2 )< 



{ 75 ,e(A-)hA) M} 



{H-f 



1 If we prefer, we can also change 9, using 9 —* e 2d 9. 



The good news is that this is (approximately) in- 
dependent of the volume. However, it is propor- 
tional to the inverse square of the mass /x (see ta- 
ble^ for a rough numerical confirmation of this). 
The probability of transmission has an exponen- 
tial dependence on d. Therefore, at small masses 
we are going to have a low acceptance rate. This 
can partially be solved by introducing multiple 
pseudo-fermion fields |1UI15) . but it still remains 
a serious issue that still needs to be resolved. We 
also expect problems changing topological sector 
at small lattice spacing. 

Zoltan Fodor and his collaborators have re- 
cently developed a novel algorithm to avoid the 
necessity of changing topological sectors 0]. It is 
easy, by continually reflecting, to keep the sim- 
ulation fixed within one topological sector. By 
starting in several different topological sectors, it 
is possible to calculate the expectation value of an 
observable with each topological sector. To cal- 
culate the total expectation value, one has to find 
the relative weighting of the various sectors. This 
can be done by measuring an observable (which is 
defined only on the topological sector wall) on ei- 
ther side of the wall. Although there are still some 
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doubts concerning the ergodicity of this method 2 , 
it does represent an interesting possibility to solve 
the problem of the topological autocorrelation at 
small masses. 

Another problem is that the fermionic force 
can be unstable if we have two small eigenvec- 
tors of opposite signs — the force is proportional 
to l/(Ai — Aa) (see equation ©). This can be 
reduced by using stout smearing JSj or an im- 
proved kernel operator to reduce the number of 
small eigenvalues [5], but these do not address the 
underlying problem. We shall discuss this issue 
further in a future publication. 

One significant advantage of overlap fermions 
is that the chiral symmetry allows us to factorise 
the squared overlap operator into the two chiral 
sectors: 

2 + 75£(Q) + e(Q)75 - 



(l + 75 )+a-(l-7 5 )+ 



1 



_(l +75 ) e (Q)(l +75 ) 



1 



-(l- 75 ) + a-(l+7 5 )- 



2 (l- 75 )e(Q)(l- 75 ) 



2 

a 



a + 1 is an arbitrary constant. Both of the fac- 
tors are positive definite, so we can use this de- 
composition to run single flavour simulations. It 
is easy to show that these two operators have the 
same non-zero eigenvalue spectrum as the overlap 
operator H (up to an unimportant sign). Zero 
modes can be included for by introducing addi- 
tional pseudo-fermion fields to generate the deter- 
minant of 1 — j-^c{Q) ■ This is only useful at large 
masses, because we have to use a polynomial or 
rational approximation to obtain the square root 
of the chiral projected overlap operator, but we 
have tested it at approximately the strange quark 
mass on small lattices (using two single flavour 
simulations) and the results for the plaquette and 
topological susceptibility agree with the 2-flavour 
HMC. This method can be used to simulate a 



2+1 (or 2+1+1 etc.) flavour theory with little 
additional effort. 

An important question is how well this algo- 
rithm scales with the volume. It has been sug- 
gested that the algorithm scales as the square of 
the volume V [§]. The work needed to perform a 
correction step is is proportional to the volume, 
and the density of small modes of the Wilson op- 
erator is also proportional to the volume. If the 
number of crossings were proportional to the den- 
sity of small eigenvalues, then this would lead to 
an 0(V 2 ) algorithm. However, it is by no means 
certain that this is the case because small eigen- 
values with opposite signs repel during the molec- 
ular dynamics. Our numerical experience is that 
although the number of crossings increases as we 
increase the volume, it scales considerably better 
than 0(V). On small lattices we observed a V 1 ' 5 
scaling for the entire HMC algorithm, although 
this needs to be checked on larger volumes. More 
work needs to be done on this area to fully answer 
this important question. 

5. RESULTS 

In figure^ we compare the small eigenvalues of 
the squared overlap operators for quenched and 
dynamical ensembles on 12 4 lattices at (approxi- 
mately) the same lattice spacing. The dynamical 
configurations were generated at a mass /i = 0.1, 
although both sets of eigenvalues were plotted 
with the massless overlap operator. It can be 
seen that fermion determinant significantly sup- 
presses the small eigenvalues of the Dirac opera- 
tor. This is, of course, to be expected, because 
small eigenvalues would reduce the determinant. 
However, it is significant because we know from 
the Banks-Casher relation that the small eigen- 
values lead to the spontaneous breaking of chiral 
symmetry. However, although the configurations 
shown in figure ^ do not posses small eigenval- 
ues, we have seen them during the MD. 3 Figure 
121 plots the eigenvalues of the overlap operator 
during one molecular dynamics trajectory. There 



2 But see the discussion in section 3 of HI . 



3 We discontinued the run which we used to generate fig- 
ure shortly after the plot was made because the lattice 
spacing was too large. We do not yet have enough data 
on our current runs to generate an improved plot. 
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Table 1 

The dependence of the potential wall on lattice size and mass [i. The ensembles were generated using 
the Luscher-Weisz gauge action. 

lattice size (5 \i < d > — < d> [i 2 



4 4 

12 4 

4 4 



7.5 
7.5 
7.5 



_M 

0.2 
0.1 
0.05 



-2.21 
-6.49 
-35.92 



0.084 
0.065 
0.090 











dynamical 

quenched >c ; 
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Figure 1. The smallest non-zero eigenvalues for the overlap operator for a \i = 0.1 dynamical ensemble 
(left), and quenched ensemble (right) on 12 4 lattices with lattice spacing ~ 0.18/m. 



were five topological charge changes during this 
trajectory. The trajectory started with a topolog- 
ical charge -1, and the eigenvalues of magnitude 
~ 0.21 are the zero modes (calculated to about 
a 10% accuracy with a mass \x = 0.1). The non- 
zero modes at either end of the plot are similar 
in magnitude to those in figure and these will 
not be responsible for chiral symmetry breaking. 
There is a considerably smaller non-zero eigen- 
value between the third and fourth topological 
charge changes. A possible interpretation of this 
figure is that we create an anti-instanton on the 
second topological charge change, and an instan- 
ton on the third. The small eigenvalue is then 

4 Although I am referring to instantons and anti- 
instantons, recent quenched calculations suggest that 
topological vacuum is not in fact dominated by instan- 
tons, but by long range topological fluctuations. But this 



be generated by the mixing between the two zero 
modes. If this interpretation is correct, then this 
presents further evidence for an underlying topo- 
logical cause for chiral symmetry breaking. 

One area in which we expect dynamical over- 
lap simulations to be particularly important is 
the measurement of the topological susceptibil- 
ity, defined as \t = (Q 2 f)/V, where Qf is the in- 
dex of the overlap operator, and V is the lattice 
volume, since overlap fermions are the the only 
lattice fermions with a well defined index theo- 
rem. In figure |21 we show plot the average value 
of the squared topological charge (QV) against 
quark mass. It proved impossible to calculate the 
lattice spacing to any accuracy on these small lat- 

argumcnt holds however the zero modes are generated: I 
only use the terms instantons and anti-instantons for the 
sake of simplicity. 
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Figure 2. The eight smallest eigenvalues of the overlap operator during a molecular dynamics trajectory. 



tices, so we have no data for \t itself. The lattice 
spacing will change with the quark mass, but not 
significantly. Although we are unable to calcu- 
late the volume, we do see the expected linear 
dependence of the topological susceptibility with 
the quark mass. 

6. CONCLUSIONS 

Dynamical overlap fermions offer an exiting 
prospect for lattice QCD at small masses. How- 
ever, the HMC is considerably harder to imple- 
ment because of the discontinuity in the .Dirac 
operator The largest problem, how to deal with 
the Dirac-delta function in the fermionic force, 
has been solved. We still have a number of smaller 
issues to resolve, such as the problem of the topo- 
logical autocorrelation at small masses and due 
to mixings between small eigenvectors of oppo- 
site signs. However, progress is being made on 
these issues. Although overlap fermions are still 
rather slow, we hope to begin large-scale sim- 
ulations in the near future. We are currently 
working on lattices with sizes up to 16 3 32 and at 
masses of approximately one third of the strange 
quark mass, and larger lattices and smaller quark 
masses should be possible on the next generation 
of computers. 
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